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I. INTRODUCTION 

Computer simulations are becoming more and more important in the study of complex 
systems. Among them many of the most interesting systems are charged systems, and as 
such naturally dominated by long range interactions. This is certainly the case for almost 
all biological systems where the electrostatic interactions play a dominant rolea, but there 
is also a wealth of technological important substances such as polvelectrolytes, where the 
charged nature is one of the key ingredient for their functionality^. The electrostatic in- 
teractions can be of monopolar origin, like interactions between proteins, DNA, or charged 
membranes, but also of dipolar nature, because all biological tissues contain water, which 
is a dipolar substance. The simultaneous appearance of dipoles and monopoles is of great 
importance, for example, for protein folding simulations^. Dipolar interactions can also be 
of technological importance, for examnle in the application of so called ferrofluids, which are 
basically dispersed magnetic particlesa. However, computer simulations of long range inter- 
actions in periodic boundary conditions are notoriously difficult to handle, and possess an 
unfavorable scaling with the amount of particles involved. The most often used method to 
compute the interactions relies on the famous Ewald sumH. In the simplest implementations 
the involved computation time grows like A^^, or at best like A^^/^, if the cutoff is opti- 
mally varied with the splitting parameteiu. The use of fast Fourier transformations (FFT) 
can further reduce the scaling to basically A^ ■ log A^. There have been quite a number of 
advances in the application of these so-called particle-mesh- Ewald techniques for Coulomb 
systems over the last few yearsQllll. An important aspect of all algorithms is the tuning in 
the sense of speed at well controlled errors. All algorithms have parameters such as the real 
space cutoff Tc, the reciprocal space cutoff fcc, the splitting parameter a, and for the mesh 
based methods there are even more. To find the optimal combination and at the same time 
control the systematic error is a formidable task, and cannot efficiently be achieved by trial 
and error. For the simulations of point charges, there have been reliable error estimates 
published for the standard Ewald summationli3, and for the mesh Ewald methods PMEli^ 
and PSMIIj. For the dipolar Ewald summation, there has been just an estimate given for 
the energy in real-space in Ref. O. For molecular dynamics (MD) simulations of dipolar 
systems, however, we need to know the errors for the forces and the torques. In this article 
we give a reliable error estimate for the energy, for the forces, and for the torques, when 
computed via the standard Ewald sum. We will show the applicability of these estimates 
by comparing them to well specified systems. Moreover we will give a detailed discussion 
on the optimization of the parameters, which will lead to the most efficient parameters for 
a predefined error in each observable quantity. This can all be done prior to the actual 
simulation, ensuring thus optimal performance at optimally controlled errors. 

The article is organized as follows: In Sec. II we briefiy review the important formulas of 
the Ewald summation for dipolar systems. The theoretical estimates of the cut-off errors in 
the Ewald sum are derived in Sec. III. In Sec. IV the estimates are compared with numerical 
accuracy measurements for specified model systems, and they are found to be very precise. 
The use of the estimates in determining the optimal parameters and their application to 
an inhomogeneous system is discussed in Sec. V and VI, respectively. Finally we end with 
some conclusions in Sec. VII. 



II. THE EWALD SUMMATION 

Consider a system of N particles with a point-dipole fj.^ at their center position r^ in 
a cubic simulation box of length L. If periodic boundary conditions are applied, the total 
electrostatic energy of the box is given by 

2 2^2^ 2^ l| r^ +n |3 |r„- + np /' ^^ 

where r^j = r^— r^. The sum over n is over all simple cubic lattice points, n = [rixL, UyL, n^L) 
with rixi Uy, Hz integers. The prime indicates that the i = j term must be omitted for n = 0. 
The slowly decaying long range part of the dipolar potential renders the straightforward 
summation of Eqn.([^) impracticable. The Ewald method provides an efficient way of cal- 
culating U, which splits the problem into two absolutely and rapidly convergent parts, one 
in real space and one in reciprocal space. The details of the method are discussed in Refs. 
^,|6|, p!5|jl6| , here we only give the final expressions 

U = f/(^') + f/(^-) + f/("^'-^) + f/(«"^/), (2) 

where the real-space f/*-^'-*, the k-space (reciprocal-space) U^'^\ the self f/(*^^-^) and the surface 
jjisurf) contributions are respectively given by: 

-^ N N I 

^^'■^ = 2 5Z 5Z E {(^^ ■ t'^)B{\ r,, + n I) - [m. ■ (r,, + n)][/x,. ■ (r,, + n)]C(| r,, + n |)}, (3) 

i=l j=l nGZS 

f^^'^ = ^ E ^exp[-(7rA:/«L)2]^5^(Mrk)(Mrk)exp(27r^k.r,,/L), (4) 

keZa,k7^0 i=l i=l 

o 3 ^ 

^"""-(27^EEm.-M. (6) 

Here the sums over i and j are for the particles in the central box and 

B{r) = [erfc(ar) + {2ar/^/^i) exp{-a'^r'^)]/r^, (7) 

C{r) = [3erfc(ar) + {2ar/^/^){3 + 2ah-^) exp(-aV2)]/r^ (8) 

The prime on the third sum in Eqn. @ also denotes that the divergent terms i = j for n = 
have to be omitted, erfc(x) := 27r^^/^ J°°exp(— t^)(it is the complementary error function. 
The inverse length a is the splitting parameter of the Ewald summation which should be 
chosen so as to optimize the performance. The form Eqn.(^) given for the surface correction 
assumes that the set of the periodic replications of the simulation box tends in a spherical 
way towards an infinite cluster and that the medium outside this sphere is an uniform 



dielectric with dielectric constant e'Ej'Ej. The case of a surrounding vacuum corresponds to 
e' = 1 and the surface term vanishes for the metallic boundary conditions (e' = oo). 

In practical calculations, the infinite sums in Eqns.(^ ^ are truncated by only taking 
into account distances which are smaller than some real space cutoff r^ and wave vectors 
with a modulus smaller than some reciprocal space cutoff k^. If r^ < L/2, the sum in real 
space (Eqn.(^) reduces to the normal minimum image convention. The double sum over 
particles in U^'^^ can be replaced by a product of two single sums which is more suitable for 
numerical calculations. 

The force F, acting on particle i is obtained by differentiating the potential energy U 
with respect to r^, i.e., 

with the real-space and k-space contributions given by: 

N / ^ 

i=l nGZ3 '^ 

- [^, ■ (r^,- + n)][/x^. ■ (r^,- + n)]D{\ n^ + n Din, + n)|, (10) 

Ff^ = ;^E E ^(M.-k)(/.,-k)exp[-(7rfc/aL)2]sin(27rk-r,,/L), (11) 

i=i kez^.k^o 

where 

D{r) = [15erfc(ar) + (2ar/v^)(15 + 10a V^ + 4a V^) exp(-aV2)]/r^. (12) 

Since the self and surface energy terms [Eqns.(|^, ^] are independent of the particle positions, 
they have no contributions to the force, unlike the Ewald summation for the Coulomb 
systems where the surface term contributes. 

The torque Tj acting on particle i is related to the electrostatic field Ej at the location 
of this particle, 

T. = M.xE, = rr) + rf)+Tr^\ (13) 



with 



E, = -1^. (14) 



and thus 

N I 



^(r) 



^ 5^ |(/x, X fi^)B{\ r,, + n I) - [/x, X (r,, + n)][/x^. ■ (r,, + n)]C(| r,, + n |)}, (15) 



i=l n&^ 



1 ^ A 

''^'^ = -L^J2 E ^(A^.xk)(AX,.-k)exp[-(7rA;/aL)2]exp(27r^k-r,,/L), (16) 

j=l keZ3,k7^0 

. N 

(surf) _ 47r ^ 

^^ --(2e' + l)L3A.^'X^r (17) 

When performing computational simulations, it is important to control the accuracy 
properly. In molecular dynamics methods the accuracy is generally estimated from the root 
mean square (rms) error in the forcesit^^ 



AF 



1^S<"^-'"^N 



1 ^ 

- ^ (F, - Ff- )^ (18) 






where Fj is the force on particle i calculated by the algorithm under investigation (here the 
Ewald sum) and F?-'^^ is the exact force on that particle. In next section, we will derive 
estimates for the rms error AF caused by cutting off the Ewald summation in real-space 
and k-space. The similar estimates are also derived for the cutoff errors in the total energy 
and torques. There are no errors involved in the self and surface contributions (Eqns(|], P, 



171)), because no cutoff operations are applied to them. 



III. EXPLICIT DERIVATION OF ERROR FORMULAS 

A. Real-Space Errors 

The real-space cutoff error in the force on particle i can be written aal3 

AFf) =1 M. I E I M, I Xff. (19) 

(r) 

The idea behind this is that the error on F^ originates from the A^ — 1 interactions of 
particle i with the other dipolar particles, and each contribution should be proportional to 
the product of the two dipole moments involved. The vector xfj gives the direction and 
magnitude of the error contribution from two unit dipoles in the simulation box, depending 
on their separation and orientations. It follows from Eqn. ([To|) that Xij is given by 

^S^ = 5Z ] [^cos^(/ii, fij) + fii cosi^ifij, r) + Aj cos^(Ai, r)]rC(r) 

r,r>rc 

— rcosi9{fl^,r) cos'&{fij,r)r^D{r) > (20) 

where fi^ and //• are unit vectors along the dipole orientations on particles i and j, f is the 
unit vector along r, and the sum in (EH) runs over all the periodic images of particle j for 



I r := r.y + n |> Tc. 'i9(/ij) Aj) denotes the angle between the vectors /ij and jlj, ^{/j-i, r) the 
angle between jj,^ and f, and so forth. 

To derive the estimate of AF^^^ , we assume that the positions and dipole moment orien- 
tations of the particles separated by distances larger than Tc are distributed randomly. This 
assumption is certainly not valid for all dipolar systems, but it is reasonable as a starting 
point. For random systems, the error contributions from different particles can be assumed 
to be uncorrelatedo 
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(21) 



where the angular brackets denote the average over all particle configurations. Obviously, 
the term ( xfj ) " the mean square force error for two unit dipoles - can no longer depend 
on i and j and is thus written as x^^^"^- Using Eqns.(|19|, ^1]), it follows that 
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(22) 



where the quantity A4 is defined as 
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To obtain the configurational average value (AF), we further assume that 
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(23) 



(24) 



which can be shown true for reasonably large systems by the law of large numbers along the 
line of reasoning of Ref . |14[ Inserting Eq. (^2]) into Eqs. ([18], ^ID we end up with the relation 
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(25) 



Using Eqn.(pOD, it follows 
X^"^^~7^ / r'^dr j smOde 

-^ Jrc Jo Jo 

— f COS 9 cos '&{fJ-', r)r^D(r] 



2tt 



[f cos 'd{il, fi) + /i cos 'd{ii , r) + fi cos 6]rC{r) 

(26) 



where jx and jx denote the unit orientation vectors of two arbitrary dipoles. The z-axis 
of the spherical coordinates (r, 9, (p) is chosen to be along the orientation fi. That gives 
'(9(/x, f) = 9 and 



COS 'd{jj,', r) = cos 9 cos 'i9(/i, /*') + sin 9 sm{}{jj,, fi) cos 0. (27) 

In the integrand of Eqn.(|26|), all the triangular functions about 'd{jj.^ fi') should be replaced 
by their configuration average values. This means that the cross terms which contain odd 
powers of cos /sin ■(9 (/i, fi') will vanish due to their zero mean value. The remaining terms 
cos^ / sin^ i){jj,, u) are replaced by their mean values of 1/2. Following the asymptotic ex- 
pansion formulaa 

/ exp{-Bx'^)f{x)dx ^ exp{-BA'^)f{A)/2BA (28) 

Ja 

which is valid when B > and d[f{x)/2Bx]/dx ^ f{x) for x > A, the complementary error 
function is approximated by 

erfc(ar) ^ exp{—a'^r'^) / y^ar. (29) 

Eqn. (|29|) is introduced into ( p6|) through (^ ^. Estimating the integral again with Eqn.(p8|) 
gives the approximate expression of x*-''''^. We obtain 

X(^)^ - L-'rfa-\^-^C! + ^^Dl - f^C^D^) exp{-2a'rl), (30) 

where the terms Cc and Dc are given by 

C, = Aa^rj + Qa^l + 3, (31) 

D^ = 8aV^ + 20aV^ + 300^^ + 15. (32) 

The resulting rms expectation of the real-space cutoff error in the forces is thus 

AFW ^ M\L'a'rlN)-'/'i^C'^ + ^D^ - ^C^D,)'/' exp{-a'rl). (33) 

6 15 15 

The derivation of the expected real-space cutoff errors in the total potential energy and 
torques follows the same way. For calculating the fluctuation of the error in total energy, 
it is noted that the interaction energy between two dipoles is evenly shared between them. 
That means the sum of {{AU^^'^Y) over all particles contains each pair contribution twice 
and thus the fluctuation of the real-space cutoff error is one half of the sumllHl. Then the rms 
value of the real-space cutoff error of the total potential energy is estimated as 

A[/M ^ M^{L'a^rlN)-'/^[\B'^ + —Cj - -B,C,Y'^ expi-a^l) (34) 

4 15 6 

with 

B^ = 2a^rl + 1. (35) 

The rms error on the torques is estimated similarly to the force as 



^ 



Eqn.(|33|, |3^ , ^6l) all contain the exponential exp(— a^r^). For sufficiently low errors, avc has 
to be larger than one, for example ar^ ~ tc for an error of exp(— tt^) ^ 5 x 10^^. If only the 
highest powers of ar^ are retained, the estimates of the real-space cutoff errors in the total 
energy, forces and torques can be reduced to 

A[/M ^ AM^a\rJ15NL^y/^ exp(-aV2), (37) 

AFW ^ 8M^a\2rl/15NLy/^ exp{-a^rl), (38) 

ArM ^ 47^2^2(^^/5^^3)1/2 exp(-aV2), (39) 

where Eqn.(p7|) is a factor of \/Q/5 slightly larger than that given in Eqn.(35) of Ref. [T2 . 
The advantage of these simplified formulas is that they reflect the dependence of the rms 
errors on a and Tc more directly and thus could be used more easily in determining the 
optimal values of these parameters. 



B. Reciprocal-Space Errors 

In deriving the estimates of the reciprocal-space (k-space) cutoff errors, we further assume 
that the radial distribution function of the particles is approximately unity at all distances. 
Following Eqn. (|ll|) , the k-space cutoff error in the force acting on particle i is given by 

AFf ^ = E E J^(^'^ ■ k)(A^. • k) exp[-(vrfc/«L)2] sin(2vrk ■ r.,/L). (40) 

i=l k,fc>fec 

Note that the diagonal term {j = i) in the sum does not depend on the positions of the 
particles. It will provide a systematic contribution to the cutoff error in k-spacell3. In 
Eqn. (|40|) this contribution equals to zero, thus there is no systematic part of the error in 
the forces. The same thing happens to the error in the torques. But for the total energy 
the diagonal terms are positive and the systematic contribution plays a dominant role in the 
cutoff error. 



The off-diagonal terms in Eqn.(^OD do depend on the positions of the particles and have 
alternating signs. The statistical approach in Sec. III. A can also be used. Similar to 
Eqn.([l9|), the off-diagonal contribution to the cutoff error in AFJ is given by 

AFlS,/ =1 ^^. I E I ^. I ^ (41) 



with 



x!?= E ^cos^(Ai,k)cos^(A„k)exp[-(7rfc/aL)2]^exp(27r2k-r/L), (42) 



where r stands for Tij and sin(27rk ■ r/L) is re-written as i exp{2TTih ■ r/L) according to the 
symmetrical character of the summation over k. k is the unit vector along k. Since the 
particles are assumed to be randomly distributed over the simulation box, the fluctuation 
of AFj. jr. can also be written as 

{{AF^Ij.r) = /i.^ E E I /^. 1 1 A^™ I ^ ■ ^S*) - f^lMV. (43) 

Again x*''^'*^ is independent of i and j. Using Eqn.(^), we have 
^{k)2^f \ / e^p[^2{7rk/aLy]k^dk smOdO cos2^(A,k)cos2^(A',k)(i0. (44) 

Choosing the z-axis of the spherical coordinates (A;, 6', 0) along the (x orientation, the same 
discussion as in Eqn.(p6D gives 

X^^^^ ~ l2M'^L-'^a^klex^[-2{'KkJaLY]/lb. (45) 

The rms expectation of the k-space cutoff error in the forces is thus 

AF^^^ ^ 8nM^L-''a{2nkl/15NY/^ exp[-{nkc/aLf]. (46) 

Here the notation AFJrl is replaced directly with AF^'^'' due to the fact of no diagonal 
contribution. 

The derivation of the off-diagonal parts of the cutoff errors in the total energy and torques 
proceeds in the same way. That the sum over ((At/^ A,)^) contains each pair contribution 
twice has also been considered in the error estimate of the total energy. The results are 
given by 

AU^J}\ ^ AM^L-^a{nkJ15Ny/^ exp[-{7rkJaL)% (47) 

Ar(^) ^ AM''L-^a{7rk,/5Ny/^ exp[-{TTkJaLf]. (48) 

At^'^^ is also used directly instead of ArLj-. 

The diagonal (systematic) part of the cutoff error in the total energy can be written as 

1 ^ A 

^Z = ir7^Y. E ;^/^^os^^(A.,k)exp[-(vrA:/«L)2] 

2_ /*00 PTT /'27r 

TT 



M^ / exp[-{TTk/aLy]k^dk / sinOde / cosH{il,k)d(l) 

Jkc Jo Jo 



LWN 
~ ^M^L-^a^k,N-^/^ exp[-{nkJaL)% (49) 

where the sum is again approximated by an integral and the asymptotic expansion formula 
of Eqn.(p8D is used to get the final estimate. The total k-space cutoff error in the total 
energy is thus 

9 



Af/('=) = Af/iS, + At/5}. (50) 

Comparing Eqn. (^91) with (^Tf ), it can be seen that the systematic part of the error is a 

1 /2 

factor of ~ Lake (^ 1) larger than the statistical part. Hence the systematic contribution 
is dominant in the k-space cutoff error of the total energy. 

Assuming that the real-space and reciprocal-space contributions to the error are inde- 
pendent, the total cutoff error in Ewald summation can be written as 

Ae = \/a0('^)' + Aew2, (51) 

where stands for U, F and r. 

IV. COMPARISON OF FORMULA WITH NUMERICAL CALCULATION 

In this section, we carry out numerical experiments to check the validity of the error 
estimates derived in the previous section. 

In order to make our measurements fully reproducible, we choose the test system to 
be consistent with the one described in Appendix D of a previous publication:ty 100 par- 
ticles randomly distributed within a cubic box of length L = 10, each of them has a 
unit point-dipole at the center. The coordinates and dipole orientations of the 100 par- 
ticles are constructed by generating 500 random numbers 7^„ (n = l,---500) between 
and 1. The first 300 random numbers give the Cartesian coordinates of the particles as 
(xj = L7^3j_2,yi = L7l3i_i,Zi = LTZsi). The remaining 200 numbers distribute the dipole 
orientations uniformly over a unit sphere's surface by the following way: 

cos0i = 27^3oo+2i-l-l.O, (52) 

(pi = 27^7^300+2^, (53) 

yielding {fi^^ = sin 9i cos (pi, fi^y = sin 9i sin (pi, fl^^ = cos 9i). The function rand which can 
be found in many C libraries is used as the random number generator. Thus the positions 



of the particles are exactly the same as that discussed in Refs. |TO|,|T^. Our unit conventions 
are as follows: lengths are measured in C and dipoles in V. Hence the unit of the energy, 
force and torque are V'^/C^,V'^/C^ and V'^/C^, respectively. 

The exact energy, forces and torques are obtained by performing a direct summation in 
real space with Eqn. (|l]) and the related derivatives. The sum over n is built up in a spherical 
way up to I ricut |= AOL. The Ewald summation calculations are carried out under vacuum 
condition (e' = 1) in order to compare the results with the direct summationllBO. The 
rms errors in the forces and torques are evaluated directly with the definition of Eqn.(|l^). 
However, it is inconvenient to do this for the energy due to the existence of the constant 
contributions ([/(^^'■^) and f/(*"'"^)), so the rms error in the total energy is simply taken as 

Af/ = ((^7-f/^^^)2)'/', (54) 
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where the angular brackets denote the configuration average. If only one specific random 
system is considered, the use of Eqn.(^) may make the result on Af/ somewhat sensitive 
to the details of the generated configuration. But it is easy to see that this sensitivity can 
be diminished by taking an average over several configurations. 

In the first step, we fix the real-space cutoff to r^ = 5.0 {L/2). The k-space cutoff kc 
changes from 2 to 12. The resulting curves for the rms errors in the forces and torques are 
plotted in Fig. 1(a) and (b) together with the analytical estimates derived in Sec. III. To 
avoid the unfavorable sensitivity to configuration details, the result on the rms error of the 
total energy [Fig. 1(c)] is taken from an average over 10 random configurations which are 
generated one after the other in the same way as the model system. Fig.l(a-c) show that for 
each kc there exists an optimal value of a which gives the minimum rms errors. For smaller 
values of a, the error contribution from the real-space is dominant, while for larger values 
the k-space contribution dominates. It can be clearly seen that the analytical estimates 
accurately predict both the real-space and k-space contributions to AF, At and AU for all 
values of a. Only for very small values of k^ deviations are observed for large a. As will be 
shown in the next section, this permits an easy way to determine the optimal parameters 
for a predefined accuracy. The discrepancy in the k-space part at very low kc is to be 
expected as the replacement of the summation over k with an integral in Eqn. ( P^ turns to 
be a rather crude approximation at this time. Since the optimal values of kc for the Ewald 
method are ranging between 7 and 25t3, this discrepancy does not affect the validity of the 
analytical estimates. As accuracy on the energy may be obscured by unimportant constant 
contributions and is sensitive to fluctuations, the rms errors on the forces and torques are 
more suitable for the accuracy measurements in MD simulations. Furthermore, since Fig.l(a, 
b) show the similar behavior of AF and At as function of a and kc (for example, they give 
the same optimal value of a for each kc), we will concentrate on discussions of AF in the 
remaining part of this paper. The same kind of results, however, have been achieved for AU 
and At in all the following cases. 

As most discussions are focused on the k-space error contributions in Fig.l by changing 
kc, we flx kc in the next step and change different values for Vc in order to further investigate 
the accuracy of the analytical estimates in predicting the real-space error contributions. 
The same model system as in Fig.l (a, b) is studied, kc is flxed to 8 and r^ is taken to 
be 3.0,3.6,4.3 and 5.0, respectively. For comparison, both Eqn.(P3D and (PSf ) are used to 
estimate AF^^\ The results in Fig. 2 show that the former formula precisely gives the real- 
space contribution to AF for all values of re- As expected, Eqn.(R8]) underestimates AF*^*^^ 



at small values of a. Since Eqns. (|37| - |39D are obtained under the assumption of arc ^ 1, the 
valid range of these simplifled formulas shifts to larger a with decrease of Vc- However, if an 
accuracy as (5 < 5 x 10~^ is required in a MD simulation so that a could not be very small 
considering the practical limitation of r^ < L/2, these estimates are almost as accurate as 
the full formulas [Eqns.(|3|), (0) and (^]. 



The last step is to demonstrate that the scaling of the rms errors with the particle 
number and dipole moment distribution is correctly given by AF ex Ai'^N~^^'^. Three 
random systems which differ only in the values of A1^ and N are investigated. The flrst 
system is the same as that studied in Fig. 1-2. The second system contains 200 particles 
among which 100 have a dipole strength of V and the other 100 have 3V. The third one 
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contains 400 particles: 100 with P, 200 with 5V and 100 with 7V. Hence their values of 
{M^,N) are respectively given by (100, 100), (1000, 200) and (10000, 400), and the ratio 
of their pref actors is thus 1 : \/50 : 50. The results of AF in Fig. 3 clearly reflect this 
scaling behavior by the constant shift of the three curves with respect to each other (note 
the logarithmic scale in the vertical axis). The analytical estimates predict the rms errors 
very precisely in all the three cases. 



V. OPTIMIZATION OF PARAMETERS 

In this section, we discuss the use of the analytical formulas derived in Sec. Ill to 
determine the optimal values of a, Tc and kc by which the required accuracy could be 
satisfied and the computation time is minimized. The detailed discussions on this subject 
can also be found in Refs. |,|1^. 

The overall computation time for computing the forces with the Ewald method is ap- 
proximately given bycJ 

r = arN'^irjLf + a^Nkl, (55) 

where the primitive overheads a^ and a^ highly depend on the implementation of the code 
and need to be found by numerical experiments. As an example, we have carried out the 
time experiments on a DEC personal workstation (CPU 433MHz) using a standard Fortran 
77 compiler. In the implementation the complementary error function and its derivative 
are calculated with table lookup and the reciprocal-space summation is optimized as in 
Refs. P,^. The linked-cell method is used to deal with the short-range forces (when doing 
simulations). The primitive overheads are then found roughly to be a^ = 2.5^s and ak = 
0.7/is. 

For a required accuracy S, the parameters a, re and kc should be chosen to minimize 



T with respect to the two constraints of the error bounds [Eqns.(0) and (|46|) ], which are 
restated as 

A = M\L'aYN)-'/\^-^C! + ^Dl - i^C.D.)V2exp(-«V,^), (56) 

= 8TTM^L'^a{27rkl/15Ny^^ exp[-{nkJaLf]. (57) 

V2 

In case oi 6 < 5x 10~^, Eqn. (|38D could be used instead of (|33D so as to show the dependence 



of the accuracy on the parameters more clearly. Eqn. ( ^6]) and ( |57D provide the qualitative 



function relations of r^ and kc with a as: rc{a) ~ —Ay/ln6/a and kc{a) ~ —By/ln5a. 
Inserting them into Eqn.(p5D and differentiating it with respect to a yields a oc N^^^ and 
thus Tc oc N~^^^ and kc oc N^^^. The minimized computation time is then proportional to 
N'^/"^ with the proportionality constant depending on the accuracy. The same results can be 



found for the Coulomb Ewald method in Refs. y,|T2|,|l3|. This can be easily understood by 
comparing Eqns.(^) and (|lB|)in Sec. Ill of this paper with Eqns.(18) and (32) in Ref. ^ 
and finding the same exponential dependences of the cutoff errors on a, Vc and kc for the 
dipolar and Coulomb Ewald summations. 
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The numerical investigation of the functional dependence of the optimal parameters on 
N and 6 are performed by using the primitive overheads obtained above. For each given 
A^ and S, we at first choose different values for r^ within the inequality r^ < L/2. For each 
Tc the parameters a and k^ are calculated by solving Eg ns . (1561) and (|57D . These values are 
then introduced into Eqn. ( |55| ) to figure out the optimal value of Vc which gives the minimum 
computation time. In calculations the size of the simulation cell is fixed to a dimensionless 
length of L = 10. The range of accuracy requirement and number of particles are chosen 
to he 6 = 10~^ to 10~^ measuring in V'^/C^ and N = 10^ to A^ = 10^ which should cover 
most of the applications. The particles are supposed to have an uniform dipole moment of 
V. The results for the optimal values of the parameters and the corresponding computation 
time per particle are shown in Fig. 4(a-d), respectively. It can be clearly seen that the 
functional dependence of the parameters and the overall computation time on A^ are just 
as discussed above. Fig. 4(c) shows that when a high accuracy is required for a system with 
a small number of particles, the predicted real-space cutoff is larger than half of the box 
length and r^ = L/2 must be used. The optimal a values hardly depend on the accuracies. 
These results are very similar to that obtained for the Coulomb Ewald summation in Ref. 
T3|, except for Vc oc N~^^^ here and Vc ex N^^^ there. This is because they considered a 



system of constant density, while we choose the volume of the simulation cell to be constant. 
VI. APPLICATION TO AN INHOMOGENEOUS SYSTEM 

So far we have only used the homogeneous random systems to test the error estimates. 
This is of course not always the case in computer experiments. Many simulations will 
encounter the problem of highly nonuniform distributions of particles. In this section we use 
an anisotropic dipolar system to investigate the influence of the inhomogeneity on the rms 
errors. 

The dipolar system we studied is an electrorheological (ER) fluid consisting of spherical 
dielectric particles with uniform diameter a dispersed in a solvent. In the initial configu- 
ration 200 particles are randomly distributed in a cubic simulation box with side length of 
lOo", which gives a typical volume fraction of particles as p ~ 0.1. Upon application of an 
external electric field, the dielectric particles are polarized. Each particle has an uniform 
dipole moment along the field direction (z-axis) at its center under the 'point-dipole' approx- 
imation. The electrostatic interaction between the particles will draw them to form chain 
or cluster structures along the field direction. This structuring process of the ER system 
in the quiescent state is studied by the Brownian dynamic simulation method described in 
our previous paperta. The final configuration, which is used for the rms error calculation, is 
obtained after a long enough simulation time when there is no obvious structure evolution 
in the system. The formation of several thick chain-like structures has been observed in the 
snapshot of the system. It is also reflected clearly by the appearance of the sharp peaks at 
the positions of r = cr, ^/Sa, 2a, Via, 3a ■ ■ ■ in the radial distribution function (RDF) go{r) 
of the final configuration shown in Fig. 5, where the RDF of the initial configuration is also 
plotted for comparison. The rms force error AF and the corresponding estimates [Eqns. ( ^3]) 
and (^)] for the final configuration are given in Fig. 6. A deviation of the estimates from the 
numerical results could be observed. At smaller values of a the prediction of the estimates 
is smaller than that given by the algorithm, while at large values of a it is larger than the 
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actual one. By studying the derivation of the error estimates, this could be understood 
qualitatively by the reason that the formation of the chain structures leads to higher local 
density of the particles and thus results in higher rms errors in the real-space parljl^. Though 
this effect also exists in the k-space part, the well ordering, or in other word well spacing, of 
particles along the z-direction leads to an enhanced k-space accuracy. This could be partially 
seen from Eqn.(|ll]) where only the z-component kz involves in the calculation due to the 
z-direction orientation of the dipole moments in this case. 

From the numerical curve in Fig.6, the optimal splitting parameter is given to he a ^ 
0.72a~^ with a corresponding AF ^ 1.74 x 10~^, while the intersection point of the real 
and reciprocal space estimates occurs at a ~ 0.71a~^, which predicts an error of AF ^ 
1.90 X 10~^ by Eqn.(^6D and (|57D . If the estimated value of a is used, this would result in 
an error of AF ^ 1.88 x 10~^, which is about 8% larger than that at the calculated optimal 
value of a. If such a safety margin has already been considered at the beginning of the 
simulation, the determination of the optimal parameters discussed in the previous section is 
still a good approach. In addition, it should be noted that the structure studied in Fig.6 is an 
extreme case in the simulations of ER fluids or similar systems such as magnetorheological 
(MR) fluids and ferrofluids. When the thermal agitation is comparable with the dipolar 
interactions, or even more a shear flow is introduced, the percolating chains or columns will 
be broken into clusters. The distribution and orientation of the smaller clusters would have 
more random characters compared with that studied in this section, and thus a smaller 
discrepancy between the estimates and actual numerical results could be expected. In any 
case, if it is suspected that the development of strongly inhomogeneities in the systems may 
lead to potential failure of the presented error formulas, some simple numerical tests as 
performed above could provide valuable informations. 



VII. CONCLUSION 

The analytical formulas for the cutoff errors in the Ewald summation for dipolar systems 
have been derived in closed form. Errors in the total energy, forces and torques are all 
considered. The high quality of these estimates has been proven in different random systems, 
and thus provides an easy way to determine the optimal tuning parameters which can give 
the expected accuracy, but minimize the computation time. Based on this, the functional 
dependence of the optimal splitting parameter a, real-space cutoff r^ and reciprocal-space 
cutoff kc on the number of particles A^ are discussed qualitatively and confirmed numerically 
by the timing experiments. Although the validity of the error estimates is subject to some 
additional requirements, such as the homogeneity of the system, a consultation of these 
formulas and a priori estimate of the optimal parameters should always be a good starting 
point in using the Ewald summation for simulations of dipolar systems. 
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FIG. 1. The rms errors (solid lines) on the forces AF (a), torques At (b) and total energy 
AU (c) as a function of the splitting parameter a. The results on AF and At are obtained for 
the model system of 100 randomly distributed dipoles, while AU is calculated from an average 
over 10 random configurations which are generated in the same way as the model system. The 
real-space cutoff is taken to be re = 5.0. The dot-dashed curves are the corresponding real-space 
error estimates achieved with Eqn.(p3|, |36| ) and (p5). The dashed curves are the k-space error 
estimates obtained with Eqn.(^, |48| ) and (|50|). The units of the energy, force and torque are 
-p2/£3^-p2/£4 g^^^ V^/jO^^ respectively. 
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FIG. 2. The rms errors (solid lines) in the forces AF for the model system. The k-space cutoff 
is taken to be kc = 8. From top to bottom the real-space cutoff varies like Tc = 3.0, 3.6, 4.3, 5.0. 
The dot-dashed curves are the corresponding real-space error estimates achieved with Eqn.(|33D. 
The results obtained with the simplified formula Eqn.(^) (the dotted curves) are also shown for 
comparison. The dashed curves are the k-space error estimates obtained with Eqn.(|46|). The unit 
of the force is V^/C^. 
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FIG. 3. Test of the Ad'^N ^'^ scaling of AF. The three sohd hnes show the rms error AF for 
systems which are different in their {M'^,N) values. From top to bottom they are characterized by 
(10000, 400), (1000, 200) and (100, 100) (the last one is the model system discussed in Fig.(l-2)). 
The real and k-space cutoffs are Tc = 5.0 and kc = 8, respectively. The dot-dashed and dashed 
curves are the corresponding real-space and k-space error estimates obtained with Eqn.(|3^) and 

1), respectively. The unit of the force is V'^ /C^. 
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FIG. 4. Optimal values of the parameters a (a), kc (b) and Vc (c) as well as the corresponding 
minimized computation time T /N (d) as a function of the number of particles. 
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FIG. 5. The radial distribution function (RDF) of the initial (dashed line) and final (solid line) 
configurations of the electrorheological system studied. The sharp peaks in the RDF of the final 
configuration reflect the formation of chain-like structures along the field direction. 
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FIG. 6. The rms force error AF for the final structure of the electrorheological system studied 
(solid line). The dot-dashed line is the error estimate for the real-space contribution [Eqn.(| 
and the dashed line is that for the k-space [Eqn.(|46|)]. The unit of the force is V'^/a'^. 
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